library(tidyverse)
library(terra)
library(sf)
library(mapview)
library(raster)
library(rgeos)
library(lubridate)
library(ggplot2)
library(exactextractr)
library(patchwoRk)
library(gridExtra)Conifer Forest Type Groups: Douglas-Fir, Ponderosa Pine, Fir-Spruce-Mountain Hemlock, Lodgepole Pine
# forest type groups and key
conus_forestgroup <- raster('data/conus_forestgroup.tif')
forest_codes <- read_csv('data/forestgroupcodes.csv')
# set crs
crs = crs(conus_forestgroup)Canadian Rockies, Idaho Batholith, Middle Rockies, Columbian Mountains - Northern Rockies
# level 3 ecoregions
l3eco <- st_read('data/us_eco_l3.shp') %>%
st_transform(., crs=crs)## Reading layer `us_eco_l3' from data source
## `G:\Other computers\My Laptop\Documents\Grad School\Research\ConiferRegeneration\data\us_eco_l3.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1250 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -2356069 ymin: 272048.5 xmax: 2258225 ymax: 3172577
## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic_USGS_version# select northern rocky mountains from level3 ecoregions
eco_select <- l3eco %>%
filter(NA_L3NAME %in% c('Canadian Rockies','Columbia Mountains/Northern Rockies','Middle Rockies','Idaho Batholith'))
mapview(eco_select)Criteria:
-1988-1991
-5000+ acres of burning
-500+ acres of high-severity
-Within selected ecoregions
->25% of selected forest types
# mtbs fire perimeters
mtbs_full <- st_read('data/mtbs_perims_DD.shp') %>%
st_transform(., crs=crs)## Reading layer `mtbs_perims_DD' from data source
## `G:\Other computers\My Laptop\Documents\Grad School\Research\ConiferRegeneration\data\mtbs_perims_DD.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 29533 features and 22 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -166.1885 ymin: 17.94736 xmax: -65.33821 ymax: 70.15893
## Geodetic CRS: NAD83# filter fires by general region of interest and date range
mtbs_select <- mtbs_full %>%
mutate(state = str_sub(Event_ID,0,2),
year = year(as.Date(Ig_Date))) %>%
filter(state %in% c("WA","ID","MT","WY","SD"),
between(Ig_Date, as.Date('1988-01-1'), as.Date('1991-12-31')))
mapview(mtbs_select,zcol="year")# function to group adjoining fire polygons to ensure patches are contiguous
group_fires <- function(mtbs_year) {
# join the polygons with themselves, and remove those that do not join with any besides themselves
combined<- st_join(mtbs_year, mtbs_year, join=st_is_within_distance, dist = 180, left = TRUE,remove_self = TRUE) %>%
drop_na(Event_ID.y)%>%
dplyr::select(Event_ID.x,Event_ID.y)
if(nrow(combined)>=1){
# partition data into that that has overlap, and that that does not
overlap <- mtbs_year %>%
filter(Event_ID %in% combined$Event_ID.x)
no_overlap <- mtbs_year %>%
filter(!(Event_ID %in% combined$Event_ID.x))
print(paste0("there are ",nrow(overlap)," overlapping polygons"))
# join all overlapping features, and buffer with radius to ensure proper grouping
overlap_union <- st_union(overlap) %>%
st_buffer(190)
# break apart the joined polygons into their individual groups
groups <- st_as_sf(st_cast(overlap_union ,to='POLYGON',group_or_split=TRUE)) %>%
mutate(year = mean(mtbs_year$year),
Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
rename(geometry = x)
print(paste0("polygons formed into ",nrow(groups)," groups"))
# join back with original dataset to return to unbuffered geometry
grouped_overlap <- st_join(overlap,groups,left=TRUE)
# arrange by the new grouping
joined_overlap_groups <- grouped_overlap %>%
group_by(Fire_ID) %>%
tally()%>%
st_buffer(1) %>%
dplyr::select(Fire_ID) %>%
mutate(year = mean(mtbs_year$year))
# add new ID to the freestanding polygons
no_overlap_groups <- no_overlap %>%
mutate(Fire_ID = str_c("Fire_",nrow(groups)+c(1:nrow(no_overlap)),"_",year)) %>%
dplyr::select(Fire_ID,year)
# join the new grouped overlap and the polygons without overlap
fires_export <- rbind(joined_overlap_groups,no_overlap_groups)
return(fires_export)
} else {
print("no overlapping polygons")
fires_export <- mtbs_year %>%
mutate(Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
dplyr::select(Fire_ID,year)
return(fires_export)
}
}# group adjacent polygons within each fire year
fires_88 <- group_fires(mtbs_select %>% filter(year == 1988))## [1] "there are 22 overlapping polygons"
## [1] "polygons formed into 7 groups"fires_89 <- group_fires(mtbs_select %>% filter(year == 1989))## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"fires_90 <- group_fires(mtbs_select %>% filter(year == 1990))## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"fires_91 <- group_fires(mtbs_select %>% filter(year == 1991))## [1] "no overlapping polygons"# join each fire year, filter by area
mtbs_grouped <- rbind(fires_88,fires_89,fires_90,fires_91)%>%
mutate(aream2 = as.numeric(st_area(geometry)),
areaha = aream2/10000,
areaac = areaha*2.471)%>%
filter(areaac >5000)# assign proportions of forest type to each fire polygon
fires_join <- st_join(mtbs_grouped,eco_select,join=st_intersects,left=FALSE,largest=TRUE) %>%
left_join(., exact_extract(conus_forestgroup,mtbs_grouped, append_cols = TRUE, max_cells_in_memory = 3e+08,
fun = function(value, coverage_fraction) {
data.frame(value = value,
frac = coverage_fraction / sum(coverage_fraction)) %>%
group_by(value) %>%
summarize(freq = sum(frac), .groups = 'drop') %>%
pivot_wider(names_from = 'value',
names_prefix = 'freq_',
values_from = 'freq')}) %>%
mutate(across(starts_with('freq'), replace_na, 0)))
# remove unnecessary columns, cleanup names
# filter to remove majority other forest types, majority not-forested cover
fires <- fires_join %>%
dplyr::select("Fire_ID","year","areaha","areaac","US_L3NAME","freq_0","freq_200","freq_220","freq_260","freq_280") %>%
rename("fire_name"="Fire_ID",
"fire_acres"="areaac",
"ecoregion" = "US_L3NAME",
"freq_df"="freq_200",
"freq_pp"="freq_220",
"freq_fs"="freq_260",
"freq_lpp"="freq_280") %>%
mutate(freq_allother = 1-(freq_0 + freq_df+freq_pp+freq_fs+freq_lpp),
freq_forested = 1- freq_0,
freq_ideal = freq_df+freq_fs+freq_lpp)%>%
mutate(across(starts_with('freq'), round,2))%>%
filter(freq_ideal > 0.25)# import all mtbs rasters via a list
rastlist <- list.files(path = "data/MTBSmosaic", pattern='.tif', all.files=TRUE, full.names=TRUE)
allrasters <- lapply(rastlist, raster)
names(allrasters) <- str_c("y", str_sub(rastlist,28,31))
# create empty dataframe
severity_list <- list()
# loop through mtbs mosasics for 1986-1991
# extract burn severity raster for all selected fires
# calculate burn severity percentages for each fire
for (i in names(allrasters)){
mtbs_year <- allrasters[[i]]
fire_year <- filter(fires, year==str_sub(i,2,5))
raster_extract <- exact_extract(mtbs_year,fire_year, max_cells_in_memory = 3e+09,coverage_area=TRUE)
names(raster_extract) <- fire_year$fire_name
output_select <- bind_rows(raster_extract, .id = "fire_name")%>%
group_by(fire_name , value) %>%
summarize(total_area = sum(coverage_area)) %>%
group_by(fire_name) %>%
mutate(proportion = total_area/sum(total_area))%>%
dplyr::select("fire_name","value","proportion") %>%
spread(.,key="value",value = "proportion")
severity_list[[i]] <- output_select
}
# combine extracted raster datasets
severity_df <- do.call(rbind, severity_list)
# join burn severity % to fires polygons
# fix naming
# filter dataset for 500 acres high severity
fires_severity <- left_join(fires,severity_df,by="fire_name")%>%
rename(noburn= "1",lowsev = "2", medsev = "3", highsev = "4",regrowth = "5", error = "6") %>%
dplyr::select(- "NaN",-"regrowth",-"error") %>%
mutate(highsev_acres = fire_acres*highsev)%>%
filter(highsev_acres > 500)
fires_select <- fires_severity %>%
left_join(.,exact_extract(conus_forestgroup,fires_severity, 'mode', append_cols = TRUE, max_cells_in_memory = 3e+08))
fires_select$mode <- as.factor(fires_select$mode)
fires_select <- fires_select %>%
mutate(fire_foresttype = case_when(mode==200 ~ "Douglas-Fir",
mode==220 ~ "Ponderosa",
mode==260 ~ "Fir-Spruce",
mode==280 ~ "Lodegepole Pine",
TRUE ~ "Other")) # plot
mapview(fires_select, zcol = "year")